Draft version February 2, 2008 

Preprint typeset using L^T^ style emulateapj v. 1 1/12/01 



cn 
O 

o 

(N 

< 

(N 
(N 

> 
oo 

o 

o 
m 
o 

i!^ ■ 

Oh' 
I 

o 



PARTICLE TRANSPORT IN TANGLED MAGNETIC FIELDS AND FERMI ACCELERATION AT 

RELATIVISTIC SHOCKS 

Martin Lemoine' 

GReCO / Institut d'Astrophysique de Paris, CNRS, 
98 bis boulevard Arago, F-75014 Paris, France 

Guy Pelletier^ 

Laboratoire d'Astrophysique de Grenoble LAOG, 
CNRS, Universite Josepti Fourier, 
BP 53, F-38041 Grenoble, France, 

Institut Universitaire de France 

Draft version February 2, 2008 

ABSTRACT 

This paper presents a new method of Monte-Carlo simulations of test particle Fermi acceleration at relativistic 
shocks. The particle trajectories in tangled magnetic fields are integrated out exactly from entry to exit through 
the shock, and the conditional probability of return as a function of ingress and egress pitch angles is constructed 
by Monte-Carlo iteration. These upstream and downstream probability laws are then used in conjunction with 
the energy gain formula at shock crossing to reproduce Fermi acceleration. For pure Kolmogorov magnetic 
turbulence upstream and downstream, the spectral index is found to evolve smoothly from s — 2.09 ± 0.02 for 
mildly relativistic shocks with Lorentz factor Fj = 2 to i ~ 2.26 ± 0.04 in the ultra-relativistic limit Fj ^ L The 
energy gain is ^ F^ at first shock crossing, and ~ 2 in all subsequent cycles as anticipated by Gallant & Achterberg 
(1999). The acceleration timescale is found to be as short as a fraction of Larmor time when Fj ^ L 

Subject headings: shock waves - acceleration of particles - cosmic rays 



1. INTRODUCTION 

Fermi acceleration at relativistic shocks is an important topic 
for understanding the formation of spectra of ultrarelativistic 
particles and radiation in relativistic flows such as those ob- 
served in active nuclei, microquasars, 7— ray bursts and pulsar 
wind nebulae (see Kirk & Duffy 2001 and references therein). 
Of particular interest is the acceleration timescale that can be 
as short as a Larmor time for relativistic Fermi acceleration; 
the smaller return probability to the shock for downstream par- 
ticles, as compared to the non-relativistic regime, is compen- 
sated by the much larger energy gain at each cycle. However 
the study of Fermi acceleration in the relativistic limit is more 
involved than in the non-relativistic regime due to the increased 
importance of the anisotropy of the distribution function (Gal- 
lant & Achterberg 1999). 

Various methods have been used to study the relativistic 
regime of Fermi acceleration (see Kirk & Duffy 2001 and ref- 
erences therein), starting with analytical estimates by Peacock 
(1981), followed by semi-analytical methods (Kirk & Schnei- 
der 1987; Gallant & Achterberg 1999; Kirket al. 2000; Achter- 
berg et al. 2001) and numerical Monte-Carlo techniques (Bal- 
lard & Heavens 1992; Ostrowski 1993; Bednarz & Ostrowski 
1998), which in spite of their differences converge to an asymp- 
totic spectral index s « 2.2 — 2.3 in the ultra-relativistic limit. 

In the present Letter, we propose a new numerical Monte- 
Carlo method of simulation of the acceleration of test particles 
at relativistic shocks. The trajectories of particles in the up- 
stream and downstream inhomogeneous magnetic fields are in- 
tegrated out exactly from the entry of each particle through the 
shock until its first return to the shock. The law of probability 
of return to the shock as a function of ingress and egress pitch 

' email: lemoine@iap.fr 

^ email: guy.pelletier@obs.ujf-grenoble.fr 

^ Unless otherwise noted, all quantities are calculated in the upstream (lab) frame; wherever needed, quantities relative to a given frame but calculated in an other will 
be marked with the subscript |, e.g., fS^^^ refers to the shock velocity measured in the downstream rest frame. 
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angles is then constructed by Monte-Carlo iteration. Finally we 
combine these probability laws, one defined for upstream and 
the other for downstream, with the energy gain formula at shock 
crossing to simulate the acceleration process. This latter use of 
the angular probability laws is similar to the method of Gallant 
et al. (2000) with some differences to be discussed below. 

The present method appears efficient and potentially more 
powerful when compared to direct Monte-Carlo simulations 
which follow each particle through its repeated shock crossings 
(e.g. Ballard & Heavens 1992; Ostrowski 1993). It allows one 
to simulate relativistic Fermi acceleration in any magnetic con- 
figuration, albeit for test particles only. We describe the method 
and numerical simulations in Section 2 and then present the re- 
sults for a planar shock with fully turbulent magnetic field both 
upstream and downstream in Section 3. 



2. NUMERICAL SIMULATIONS 

The hydrodynamic jump conditions at an adiabatic strong 
shock, neglecting magnetic fields, are given in Blandford & 
Mc Kee (1977), Kirk & Duffy (2001) and Gallant (2002). The 
shock Lorentz factor is Fj (upstream or lab frame^), and the 
relative Lorentz factor Fj between upstream and downstream 

Id 



rs|d(l ~ PsPs\d)- The downstream Lorentz factor Fj 



(and Fi) can be obtained as a function of Fs from the rela- 
tions derived from the shock jump conditions given in Gallant 
(2002). In particular, in the ultra-relativistic limit Fs +00, 
one finds the well-known results /3s|d ^ 1/3 and Fi ^ Fj/v^. 

We conduct our simulations in two steps. We first perform 
Monte-Carlo simulations of particle propagation in a given 
magnetic field structure following Casse, Lemoine & Pelletier 
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(2001) (wherein one may find more details on the numerical 
procedure). These simulations are carried out separately in the 
downstream and in the upstream rest frames. It is possible to 
set up any magnetic field structure including regular and tangled 
components, but in the following, for the sake of simplicity, we 
restrict ourselves to the case of pure Kolmogorov turbulence 
in both downstream and upstream rest frames. The equations 
of motion of each particle are integrated out exactly, the mag- 
netic field being calculated at each point of the trajectory as a 
sum of plane waves, using 200 modes spaced logarithmically 
on three decades of wavelength and whose wavevector direc- 
tions are drawn at random^. Ostrowski (1993) used a similar 
technique to construct the magnetic field albeit with 3 modes 
only, while Ballard & Heavens (1992) used three-dimensional 
FFT methods (see Casse, Lemoine & Pelletier 2001 for a com- 
parison of these methods). 

The trajectories are integrated over 100 scattering times up- 
stream and 1000 scattering times downstream in order not to 
miss possible late returns to the shocks. The laws of return 
probability as a function of ingress and egress pitch angles are 
then constructed in the following way. We draw at random a 
point along a simulated trajectory which defines the point of en- 
try through the shock. The ingress pitch angle cosine to shock 
normal /x' is recorded, the trajectory scanned to find the point 
of exit through the shock, and the egress pitch angle cosine jjf 
is then recorded. Iteration of the above then yields the desired 
law of conditional return probability 7^(/i'; ^f) which gives the 
probability for a particle entering with a pitch angle cosine /i' 
to return to the shock with a pitch angle cosine jjp. The simu- 
lations also give a direct measurement of the return time to the 
shock as a function of pitch angles. 

Once the upstream and downstream laws of return prob- 
ability, respectively 'Pu(Mu;Mu) ^d(Md'A*d) known, 
the simulation of the acceleration process itself can be per- 
formed as follows. We denote by /^"^' (/Xd, ed) the distribution 
function of particles that enter the shock to downstream and 
that have experienced 2n+l shock crossings, and /u"(/^u,eu) 
similarly upstream: particles are upstream for an even num- 
ber of shock crossings and represents the injection popu- 
lation. These distribution functions are normalized to the to- 
tal number of particles injected such that, in the absence 
of escape from the acceleration site, after 2n shock cross- 
ings N = J d/iudeu/u"(/Ltu) eu)> and after 2n+l shock crossings 

= /d/idded/j"''''(/id,ed)- The conservation of particle num- 
ber at shock crossing u ^ d and Lorentz transforms of pitch 
angles and energies imply: 



/d"+'(ed,Md)dMdded = 



Cd = r,(l-A/iS)eu, 



and one obtains a similar system for shock crossing d ^ u: 



d/iydeu 

(1) 
(2) 



/u (eu,Mu)dMudeu = 



M^)/r-He-d;/i'd) 



Mi. = 4^' e„ = r,(l+/iM^)ed, 



d/ijded 
(3) 
(4) 



Using a higher number of modes does not modify the results shown here as we have 



The terms within brackets in Eqs. (1) and (3) correspond to 
the distribution function upon exit from upstream and down- 
stream respectively. The return probabihty to the shock /ret(/Ud) 
as a function of ingress pitch angle can be obtained as: 
^ret(Md) = /d/id^d(/Xd;Md)- The corresponding return proba- 
bility for upstream is obviously unity. After each cycle u 

J ^ M, a fraction f^:+'{e) = /d/xj,[l - P,et(Mli)]/d"+'(M!i;e) 
of the particle population has escaped downstream and ac- 
cumulates to form the outgoing particle population /out(e) = 
X^"lJ^°°/out^^(e). By following each shock crossing, and us- 
ing Eqs. (1),(2),(3), (4) one can follow the evolution of /d, /u 
and /out, starting from a mono-energetic and isotropic initial 
injection distribution upstream. The distribution /out(e) even- 
tually provides the accelerated particle population. A similar 
formal development of the acceleration process by repeated 
shock crossings has been proposed independently by Vietri 
(2002): the flux of particles crossing the shock in the station- 
ary regime, noted Ji„ in Vietri (2002) is related to the above as 

•^in = Cjy^Zo'^ /d"^' with C a normahzation constant. 

This method assumes that the angular probability laws do 
not depend on rigidity. This is true in the diffusive limit but one 
might expect some weak dependence to appear in the relativis- 
tic limit Fs » 1 . Indeed we have found numerically such a weak 
dependence of and Vu on the particle rigidity. However it 
remains weak, and the change in V averages to a few percent 
when the rigidity changes by two orders of magnitude. In terms 
of energy spectral index s, this dependence introduces an error 
of Ss = ±0.02 for = 2 up to Ss = ±0.04 for = 100. There- 
fore in the following we neglect the dependence on rigidity but 
keep the above errors as uncertainties on our results. Note that 
one can in principle incorporate this dependence on rigidity in 
our method at the expense of having to calculate the probabihty 
laws V for a wide range of values of the rigidity. 

The present technique has significant advantages when com- 
pared to standard Monte-Carlo techniques which follow the 
particle trajectories on both sides of the shock through the 
whole acceleration process. Indeed, provided one neglects the 
dependence on rigidity of V, one can simulate the trajectories 
of particles of high rigidity only (near the end of the resonance 
range) which are must faster to integrate than the trajectories 
of particles of small rigidity since the ratio of scattering time 
to Larmor time decreases with increasing rigidity. The direct 
Monte-Carlo methods also suffer from the problem of a small 
dynamic range of the magnetic fields as compared to the wide 
dynamic range of particles momenta. The present method also 
offers a significant gain in signal as will be obvious in the fol- 
lowing. Finally our method differs from Gallant et al. (2000) 
as they use Ito differential techniques to simulate scattering 
downstream and analytical methods for scattering in a regu- 
lar magnetic field upstream assuming Fg 1. In contrast, the 
present simulations can be applied to any shock Lorentz factor 
and magnetic field configuration. Furthermore they use Monte- 
Carlo methods to simulate the acceleration process after con- 
structing the laws of return probability while we directly fold 
over repeatedly the probabihty distributions in conjunction with 
the energy gain formula. 



3. RESULTS 

The downstream return probability to the shock as a function 
of ingress pitch angle cosine is shown in Fig. 1 for Fj = 2, 100. 

checked. 
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The return probability for Fs = 100 is an asymptote which is 
reached to within a percent as early as Fs = 5. In Fig. 2, 
we show the average energy gain {g) = (e/)/(e,) per cycle 
u^d^u (diamonds) and d^u^d (triangles) for Fs = 100. 
This energy gain is defined as the ratio of the average ener- 
gies at the end (e/) and at the beginning of the cycle (e,), with 
(e) = / d/idee/(/i,e)/ / d;ude/(/i,e). The average energy gain 
in each cycle subsequent to first shock crossing is ~ 1 .93 per 
cycle u-d-u for Fs = 100, and this asymptotic value is reached 
immediately after the first cycle. This is a rather dramatic con- 
firmation of the analytical expectations of Gallant & Achterberg 
(1999) and Achterberg et al. (2001) which had argued that only 
the first cycle should yield a gain w F^ since the anisotropy 
of the distribution function upstream is so pronounced that the 
gain in subsequent cycles is reduced to « 2. 

An example of the spectrum of accelerated particles for Fs = 
100 after 20 cycles m — > J — > m is shown in Fig. 3; the thin Unes 
in this figure show the fractions of particles fl"^^ that escape 
after 2n + \ shock crossings. One clearly sees in this figure the 
piUng up of populations of particles of ever decreasing size (due 
to finite escape probability) and ever increasing energy which 
gives rise to the accelerated population /out — Sn/out^'- The 
spectral index of the escaping population for Fs = 100 is here 
i = 2.26 ± 0.04 (incorporating the error due to dependence of 
V on rigidity), in excellent agreement with previous results by 
Bednarz & Ostrowksi (1998), Kirk et al. (2000) and Achterberg 
et al. (2001). 

Finally, in Fig. 4, we give the average return probabil- 
ities (open squares), average asymptotic energy gains (dia- 
monds) and spectral indices (filled circles) for values of Fs 
comprised between 2 and 100. The average return prob- 
abilities shown in this figure are the return probabilities 
shown in Fig. 1 weighted by the corresponding asymptotic 
downstream ingress pitch angle distribution, i.e. {Pkx) = 

lim„^+co/'Pret(/i!j)/d(Md)d/id///d(Md)dMd- A naive un- 
weighted average of the return probability shown in Fig. 1 for 
Fs = 100 would give 0.33, whereas the weighted average gives 
0.40: the difference is directly related to the strong anisotropy 
at shock crossing. 
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Fig. 1 . — Probability of return to shock Pmt(m) downstream as a function 
of ingress downstream pitch angle cosine for Fs = 2, 100. Note that the 
probabihty is defined for — 1 < /ij < Aid due to shock crossing conditions on 
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Fig. 2. — Averaged energy gain per cycle m — > — » m (diamonds) and 
d^u—*d (triangles) for Fs = 100. 
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Fig. 3. — Spectrum of particles escaping downstream (thick line) as a func- 
tion of momentum after 20 cycles for Fs = 100; in thin lines, the spectra of 
particles escaping downstream after each cycle. 
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Fig. 4. — Average return probabilites (open squares), average energy gain 
per cycle (dopen diamonds) and spectral slope (filled circles) as a function of 
Fs. The dotted line shows the approximation to the spectral slope given by the 
Bell formula using the average return probabilities and energy gains (see text). 

The standard non-relativistic formula for the (energy) spec- 
tral index s (Bell 1978), s = 1 - log ((Pret» /log ((g)) , with {g) 
the average energy gain, is in relatively good agreement with 
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the slopes obtained, provided one uses the weighted average 
for the return probabihty as above, see Fig. 4. A more general 
formula which includes relativistic effects has been proposed 
by Vietri (2002): {P,^i){g'-^) = I. One can derive this for- 
mula and variants of it by using our Eqs. (1),(2),(3), (4). For 
instance, insert Eq. (3) into Eq. (1) then sum over n (shock 
crossing number) to go to the stationary regime and consider 
an energy range where e ^ eo, eo being the maximal injection 
energy. There one expects that X^n/d"^' ^ ■'(?!)(/Xd), i.e. the 
distribution factorizes out in an energy power law times a func- 
tion of pitch angle, and indeed this property is verified numer- 
ically to high accuracy. Then one introduces the energy gain 
per cycle = ej/fd and integrates over both sides of 

the equation. Finally, dividing one side by the other yields the 
following relation, which is a variant of the formula of Vietri 
(2002): 

/ ''d/xj, / df,lPrMVn{l4,^J},)g'-\l4,^J,i,) = i . (5) 

In this equation, ■Pu(Md' Md) simply corresponds to Vn expressed 
in terms of downstream pitch angle cosines, and ^ret(Md) = 
/ d^Jj7'd(A*d'Md)- Equation (5) is indeed verified to within the 
numerical noise (< 1%). 

Our simulations also provide a direct measurement of the ac- 
celeration timescale, which can be defined as the cycle time 
in the upstream rest frame when Fs 3> 1: facc(e) ~ fu|u(e) + 
rr?d|d(e/rr), where ?u|u(e) and ?d|d(e) denote the upstream 
and downstream return times measured in their respective rest 
frames for a particle of energy e. For f^iij we find to within 

the noise of the simulations: t^^^ ~ 1 -S/^sJtJ^scattld' with tscM\d 
the scattering time downstream. The scattering time is given 
as a function of Larmor time in Casse, Lemoine & Pelletier 
(2001), and for pure turbulence one finds: 4catt/?L — OAp~^/^ 
for p < 0.1 and tscm/th ^ 2 for 0.1 < p < 1. Here p = feminn. 
denotes the rigidity, with kmin the smallest wavenumber of the 
magnetic field modes. The above result for ?d|d can be under- 
stood as follows. In the non-relativistic limit one derives f(j|(j ~ 

(2/3/32 J(l - (P,e,))/(A-et> and we find (P,et> - (1 - /3s|d)/2 in 
the relativistic limit, which gives fd|d ~ fscatt|d//3s Id- 
Gallant & Achterberg (1999) have conjectured ?u|u ~ ^L|u/rs 



for Fs 1, corresponding to deflection by an angle of or- 
der 1 /Fs in a regular magnetic field. We confirm this result 
up to a small residual dependence on scattering time/rigidity: 
tu\u ~ 5F~'p~'"^fL|u for Fs > 5. In the mildly relativistic case, 
for Fs = 2, we find « fscattlu^ which indicates that the particles 
have time to wander before being caught by the shock in this 
case. 

The final acceleration time depends on how the magnetic 
field downstream Bd is related to that upstream Su (Gallant & 
Achterberg 1999). If one assumes that i?d — ^b^u with F5 ^ F^, 
and the turbulence is isotropic downstream but the length scale 
has been contracted by a factor ~ F^, i.e. fcmin|d ~ rB^min|u> one 
finds, for p^j^^ ~ 1, face ^ (SF^" '^ +9)fL|u/rs- The acceleration 
timescale can thus be as short as a fraction of Larmor time in 
the ultra-relativistic limit. 

This is a most interesting aspect of relativistic Fermi accel- 
eration, as it implies that the particles can reach the energy 
confinement limit — ZeTBr when the acceleration is lim- 
ited by expansion losses or by the age of the shock wave. Here, 
r = mm{l,t/Kc), I is the size of the accelerating region, t the 
age or losses timescale and k = facc/^L < 1, all quantities being 
expressed in the comoving frame. This is particularly relevant 
for the generation of ultra-high energy cosmic rays in relativis- 
tic winds such as 7— ray bursts (Waxman 1995, Vietri 1995, 
Gallant & Achterberg 1999, Gialis & Pelletier 2003). In par- 
ticular our results for Fs = 2 are of direct relevance to the ac- 
celeration of protons and electrons in internal shocks of 7— ray 
bursts, while the results for Fs » 1 can be applied directly to 
the external shock model of 7— ray bursts. 

To sum up, our simulations confirm the results of Bednarz 
& Ostrowski (1998) concerning the spectral index and those of 
Achterberg et al. (2001) concerning the spectral index, the ac- 
celeration time scale and energy gains. Very recently, Ellison 
& Double (2002) found a similar index by taking into account 
the backreaction of cosmic rays on the shock structure. In the 
near future, the present method will be applied to more gen- 
eral magnetic field configurations including parallel, transverse, 
subluminal and superluminal shocks. 



We would like to thank Y. Gallant for fruitful discussions. 



REFERENCES 



Achterberg, A., Gallant, Y. A., Kirk, J. G., & Guthmann, A. W. 2001, MNRAS, 

328, 393 

Ballard, K. R., & Heavens, A. F. 1992, MNRAS, 259, 89 

Bednarz, J., & Ostrowski, M. 1998, Phys. Rev. Lett., 80, 3911. 

Bell, A. R. 1978, MNRAS, 182, 147. 

Blandford, R., & McKee, C. 1977, Phys. Fluids, 19, 1130. 

Casse, F, Lemoine, M., & Pelletier, G. 2001, Phys. Rev. D, 65, 023002. 

EUison, D., & Double, G. 2002, Astroparticle Physics, 18, 213. 

Gallant, Y. A., & Achterberg, A. 1999, MNRAS, 305, L6. 

Gallant, Y. A. et al. 2000, in Gamma-Ray Bursts, eds. R. M. Kippen et al., 

AIP Conference Series, Vol.526, p.524 (New-York: American Institute of 

Physics) 



Gallant, Y. A. 2002, to appear in Relativistic Flows in Astrophysics, eds. A. 
W. Guthmann et al.. Lecture Notes in Physics (Berlin: Springer- Verlag), 

arXiv:astro-ph/02 0124 3 
Gialis, D. & Pelletier, G., Astropart. Phys. in press, 

arXiv: astro -ph/03022 31. 
Kirk, J. G., & Schneider, R 1987, ApJ, 315, 425 
Kirk, J., & Duffy, P 1999, J. Phys. G, 25, R163. 

Kirk, J., Guthmann, A., Gallant, Y., & Achterberg, A. 2000, ApJ, 542, 235. 

Ostrowski, M. 1993, MNRAS, 264, 248 

Peacock, J. 1981, MNRAS, 196, 135. 

Vietri, M. 1995, ApJ, 453, 883. 

Vietri, M. 2002, arXiv:astro-ph/ 02 12352. 

Waxman, E. 1995, Phys. Rev. Lett., 75, 386. 



